/**
 * Autor:   Lukas Brzobohaty, xbrzob06@stud.fit.vutbr.cz
 *          David Marecek, xmarec10@stud.fit.vutbr.cz
 * Soubor:  main.cc
 * Projekt: Implementace modelu zalozeneho na metode Monte-Carlo
 * Predmet: IMS
 * Popis:   
 * Poznamka:
 **/

#include "simlib.h"
#include "generator.h"
#include <iostream>
#include <vector>
#include <math.h>
#include <time.h>

using namespace std;

#define POCET_KONSTANT 9
#define POCET_INTERVALU 10
#define POCET_GEN_VZORKU 100

/**
* Funkce pro vypocet prumeru a rozptylu
*
*
**/

void vypocet_prumer_rozptyl(int konst, double* vyberovy_prumer, double* vyberovy_rozptyl)
{
	y = rand();	// vygenerovani pocatecni hodnoty seed
	
	vector< vector<double> > pole_vzorku(POCET_INTERVALU);

	for(int i = 0; i < POCET_GEN_VZORKU; i++)	// generovani vzorku
	{	
		double vzorek = 0;
		int index = 0;
		if(konst == 0)	// b_0 - rozsah 0.03 - 0.1
		{
			double rozsah = (0.1 - 0.03) / POCET_INTERVALU;
			vzorek = Generator()*(0.1-0.03)+0.03;
			for(int i = 1; i < POCET_INTERVALU +1; i++)
			{
				if((i * rozsah + 0.03) > vzorek)
				{
					index = i - 1;	
					break;				
				}
			}
		}
		else if(konst >= 2 && konst <= 4)	// alfa_1, alfa_2, alfa_3, beta_1, beta_2, beta_3 - rozsah 0.008 - 0.01
		{
			double rozsah = (0.01 - 0.008) / POCET_INTERVALU;
			vzorek = Generator()*(0.01-0.008)+0.008;
			for(int i = 1; i < POCET_INTERVALU +1; i++)
			{
				if(i * rozsah + 0.008 > vzorek)
				{
					index = i - 1;	
					break;			
				}
			}
		}
		else	// rozsah 0 - 1
		{
			double rozsah = (1 - 0) / (double)POCET_INTERVALU;
			vzorek = Generator();
			for(int i = 1; i < POCET_INTERVALU +1; i++)
			{
				if(i * rozsah > vzorek)
				{
					index = i - 1;				
					break;
				}
			}
		}

		pole_vzorku[index].push_back(vzorek);
	}

	for(int i = 0; i < POCET_INTERVALU; i++) 
	{
		int pocet_vzorku = 0;
		double soucet_vzorku = 0;
		for(int j = 0; j < pole_vzorku[i].size(); j++) // vypocet vyberoveho prumeru
		{ // suma vzorku
			soucet_vzorku += pole_vzorku[i][j];
			pocet_vzorku++;
		}
		vyberovy_prumer[i] = soucet_vzorku / pocet_vzorku;
		//cout << vyberovy_prumer[i] << endl;
		
		double soucet_rozptylu = 0;
		
		for(int j = 0; j < pole_vzorku[i].size(); j++) // vypocet vyberoveho rozptylu
		{
			soucet_rozptylu += (pole_vzorku[i].back()-vyberovy_prumer[i])*(pole_vzorku[i].back()-vyberovy_prumer[i]);
		}
		
		vyberovy_rozptyl[i] = soucet_rozptylu / pocet_vzorku;
	}
}

double gen_konst_prumer[POCET_INTERVALU][POCET_KONSTANT];
double gen_konst_rozptyl[POCET_INTERVALU][POCET_KONSTANT];

struct kulm_prutok
{
	Integrator Qbmax;
	kulm_prutok(double M, double b_0, double H_k, double B, double H_d):
		Qbmax(M*(b_0*pow(H_k,3/2) + 0.4*((B-b_0)/H_d)*pow(H_k,5/2))) {}
};

kulm_prutok Q(1.5, gen_konst_prumer[0][0], 0.86, 0.2, 0.86);

// Output sampling
void Sample() {
    Print("%f %f\n", T.Value(), Q.Qbmax.Value());
}
Sampler S(Sample, 0.01);        // Output step 

/**
* Hlavni funkce
*
*
**/
int main()
{
	srand(time(NULL));

	/*   
  * b_0 - počáteční šířka nátrže
  * s_0 - průměrný sklon svahu natrze 
  * alfa_1 alfa_2 alfa_3 beta_1 beta_2 beta_3 - koeficienty eroze
  * Q_in - pritok vody do nádrže
  * Q_in_t - casova zavislost pritoku do nadrze
	*
	* 		b_0 s_0 alfa_1 alfa_2 alfa_3 beta_1 beta_2 beta_3 Q_in
	*    i	 0	 1	  2		 3	    4	   5	  6		 7	    8
	* j
	* 1
	* 2
	* 3
	* ...
	* ...
	* ...
	*/

	double vyberovy_prumer[POCET_INTERVALU];
	double vyberovy_rozptyl[POCET_INTERVALU];
	
	for(int i = 0; i < POCET_KONSTANT; i++) 
	{
		for(int j = 0; j < POCET_INTERVALU; j++) 
		{
			vypocet_prumer_rozptyl(i, vyberovy_prumer, vyberovy_rozptyl);
			gen_konst_prumer[j][i] = vyberovy_prumer[j];
			gen_konst_rozptyl[j][i] = vyberovy_rozptyl[j];
		}
	}
	
	// vypis vygenerovanych hodnot
/*	Print("b_0\t\t s_0\t\t alfa_1\t alfa_2\t alfa_3\t beta_1\t beta_2\t beta_3\t Q_in\t Q_in_t\n");
	for(int i = 0; i < POCET_INTERVALU; i++)	
	{
		for(int j = 0; j < POCET_KONSTANT; j++)
		{
			Print("%f\t ",gen_konst_prumer[i][j]);
		}
		Print("\n");
	}*/

	double H_0 = 0.84;
	double Z_0 = 0.81;
	double M = 1.5;
	double M_t = 1.35;
	double A_s = 60.932;
	double b_0 = gen_konst_prumer[0][0];
	double s_0 = gen_konst_prumer[0][1];
	double Q_in = 0.0232; //gen_konst_prumer[0][8];
	double Q_f = 0;
	double alfa_1 = gen_konst_prumer[0][2];
	double alfa_2 = gen_konst_prumer[0][3];
	double alfa_3 = gen_konst_prumer[0][4];
	double beta_1 = gen_konst_prumer[0][5];
	double beta_2 = gen_konst_prumer[0][6];
	double beta_3 = gen_konst_prumer[0][7];
	double delta_t = 0.01;

/*	for(int i = 0; i < POCET_INTERVALU; i++)
	{
	}*/

	double Q_b = M * b_0 * pow(H_0 - Z_0,3/2) + M_t * s_0 * pow(H_0 - Z_0, 5/2);
	H_0 = (delta_t / (A_s * H_0)) * (Q_in - Q_b - Q_f) + H_0;

	double A_b = b_0 * (H_0 - Z_0) + s_0 * pow(H_0 - Z_0, 2);
	double v_0 = Q_b / A_b;

	Z_0 = - delta_t * alfa_1 * pow(v_0, beta_1) + Z_0;
	
	if(Z_0 > 0)
	{
		b_0 = delta_t * alfa_2 * pow(v_0, beta_2) + b_0;
	}
	else if(Z_0)
	{
		b_0 = delta_t * alfa_3 * pow(v_0, beta_3) + b_0;
	}

	cout << Q_b << endl;

//	Print("b_0_min: %f\t b_0_max: %f\n", gen_konst_prumer[0][0], gen_konst_prumer[29][0]);
//	Print("alfa_1_min: %f\t alfa_1_max: %f\n", gen_konst_prumer[0][2], gen_konst_prumer[29][2]);
//	Print("beta_1_min: %f\t beta_1_max: %f\n", gen_konst_prumer[0][5], gen_konst_prumer[29][5]);

	/*SetOutput("ims.dat");
  	Print("# Time   Qbmax \n");
	Init(0, 200);
	Run();
	SIMLIB_statistics.Output(); // print run statistics
*/
	//cout << endl << "Byl vygenerovan soubor pseudonahodnych cisel o velikosti " << pocet_vzorku << " cisel" << endl;
	//cout << "Rozdeleno na " << POCET_INTERVALU << " intervalu" << endl;

	// Vypocet 
	/*
	while(pocet_iteraci > 0)
	{
		vemu pseudonahodne cisla
		dosadim do vzorce
		vysledek
		ulozim
		
	}	
	*/
	//cout << "Vysledne grafy" << endl;
}
